function [ ] = calculateB0FieldMap_3t(b0_te_1, b0_te_2, b0_map, span, delta_t)
%CALCULATEB0FIELDMAP_3T Script that produces a B0 map from 3T data (MINC)
%   calculateB0FieldMap_3t is a function that produces a B0 map from 3T data
%   from data formated with the MINC. It calculated the differences of two
%   phase images from different TEs, unwraps the difference, then scales it
%   to a frequency B0 map.
%
%   *Input Parameters*
%   
%   b0_te_1 : Filename (string) for the first echo. Must be located in the 
%   current directory.
%
%   b0_te_2 : Filename (string) for the second echo. Must be located in the 
%   current directory.
%
%   b0_map : Filename (string) of the output b0 map. Will be created in the
%   current directory.
%
%   span : Range of pixel values (Maximum pixel value - minimum pixel
%   value). Default value is 8192.
%
%   delta_t : Difference in echo times between second and first echo, in
%   seconds.
%
%
%   Author: Mathieu Boudreau
%   Date Created: March 20th 2013

%% Function argument checks
%
if (nargin < 2)
    fprintf('Provide TWO phase maps\n')
    return
end

if (nargin < 3) 
    fprintf('Provide name of output file.\n')
    return
end


if (nargin < 5 )
	m_span = 8192;
    m_delta_t =4.48e-3;
    fprintf('\nSpan and delta_t values were not in the function arguments.\n\nSpan and delta_t were set to default values\n');
    fprintf('\nm_span = %d \nm_delta_t = %d\n', m_span, m_delta_t);
end

if (nargin == 5)
    m_span = span;
	m_delta_t = delta_t;
end

%% Preliminary calculations / check
% 

if (m_span == 8192) 
    m_scale = 1/(4096*m_delta_t); % Converts pixel value to delta frequency.
    m_max = 4096;    
else
    fprintf('\nSpan not supported: %d\n', m_span)
end

%% Check existence of phase image files and b0 directory
%

if (~exist(b0_te_1,'file'))
    fprintf('\nPhase image 1 not found: b0_te_1\n')
    return
elseif (~exist(b0_te_2,'file'))
    fprintf('\nPhase image 2 not found: b0_te_2\n')
    return
else 
    fprintf('\nAll required maps found.\n')
end

%% Process phase maps
%

% Calculate phase difference
if (exist(sprintf('b0_map.mnc'),'file'))
    sound([sin(5*[0:0.1:500]).*sin(4*[0:0.1:500]),sin(7*[500:0.1:1000]).*sin(6*[500:0.1:1000]),sin(9*[1000:0.1:1500]).*sin(8*[1000:0.1:1500])]);    
    fprintf('\nb0_map.mnc already exists. Press "control-X" to stop.\n')
    fprintf('\nFiles will be overwritten in 5 seconds.\n')
    start = tic;
    while (toc(start)<10)
      pause(0.001)
    end    
end

system(sprintf('mincmath -sub -clobber %s %s b0_phase_diff_mri.mnc', b0_te_2, b0_te_1))

fprintf('\nPhase difference calculated.\n')

% Un-wrap phase difference
fprintf('\nUn-wrapping phase difference map.\n')

system(sprintf('minccalc -clobber -expression ''A[0]>''%d''?A[0]-''%d'':A[0]'' b0_phase_diff_mri.mnc b0_phase_diff_unwrap_tmp_mri.mnc', m_max, m_span))
system(sprintf('minccalc -clobber -expression ''A[0]<-''%d''?A[0]+''%d'':A[0]'' b0_phase_diff_unwrap_tmp_mri.mnc b0_phase_diff_unwrap_mri.mnc', m_max, m_span))

% Output B0 map
fprintf('\nCalculating B0 map...\n')
system(sprintf('mincmath -clobber -scale -const %s b0_phase_diff_unwrap_mri.mnc %s', m_scale, b0_map))
fprintf('\nB0 map calculated.\n')

end

